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INTRODUCTION 

SOILSIM is a digital model of energy and moisture fluxes in the soil and above the soil surface. 
It simulates the time evolution of soil temperature and moisture, temperature of the soil surface 
and plant canopy the above surface, and the fluxes of sensible and latent heat into the atmosphere 
in response to surface weather conditions. The model is driven by simple weather observations 
including wind speed, air temperature, air humidity, and incident radiation. The model is intended 
to be useful in conjunction with remotely sensed information of the land surface state, such as 
surface brightness temperature and soil moisture, for computing wide area evapotranspiration. 
The model was developed by Peter Camillo and reported on in a number of papers (Camillo and 
Schmugge, 1983; Camillo, Gurney and Schmugge, 1983; Camillo, 1986;). Prior to his death in 
1988, Camillo had substantially modified the program to include an option for a ’’two layer" 
model structure which allowed separate values to be solved for the soil surface and canopy 
temperatures. A working version had been developed of the revised model which was used with 
data from the HAPEX exercise (Camillo, 1988). However, this FORTRAN code was not well 
described and was not usable by the uninitiated. An early version of the model was carefully 
documented by Camillo and Schmugge (1981). Some of the substantial changes made to the 
FORTRAN code after this time were partially described informally by Camillo in several draft 
documents. This description has been prepared from Camillo’s original and revised 
documentation to the extent possible. Description of the scientific and mathematical basis of the 
modifications and additions to the FORTRAN code had to be worked out by reverse engineering 
the code in conjunction with the references with which he was working. In the present 
documentation, Camillo’s original text is retained where it applies. The section describing the 
canopy modelling was written by Field. Procedures have been added so that simulator now 
computes local mean time in addition to simulator time. This was added to rectify the times of 
simulator output and observations. Also the declination of the sun is now computed from data 
on latitude and day of the year. The list of references has been updated to include additions to 
the model and results from model simulations. 

In brief, the approach integrates a pair of differential equations for soil heat and moisture flux 
in each soil layer to satisfy certain deep soil and soil surface boundary conditions. The chief 
boundary condition at the soil surface is the total flux of heat between the soil and air above the 
soil. The fluxes of radiation, sensible and latent heat above the soil surface, whether bare or 
covered with a plant canopy, are computed iteratively as responses to the forcing induced by the 
observed weather observations. The solved above-surface "boundary fluxes" are such that their 
sum at the soil surface equals the surface value of the solved soil heat flux. The soil surface 
temperature must satisfy both the below-surface and above-surface fluxes. When a canopy is 
present, temperatures for plant leaves and canopy air are calculated. These may be either the 
same as the surface temperature (a one layer model) or different from it (a two layer model). 
Comparing surface temperature or canopy temperature with observations provides a test of the 
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quality of energy budget solution. Alternatively, one may compare observations of soil moisture 
with calculated values. 

The model was originally written to run on an IBM 3060. The model has been ported to the DEC 
VMS environment to run on VAX computers by Karen Humes, University of Arizona. It was 
ported from this version to the UNIX environment on a SUN 386i workstation by John Kuhn and 
Richard Field at the University of Delaware Center for Remote Sensing. Rather extensive 
modifications were necessary were to provide file interfacing to UNIX and the windowing 
environment of the SUN operating system (Sun OS 4.0). The program has been satisfactorily 
tested on both the Sun 3 and Sun 386i computers. In what follows only those changes to 
variables which appear in the FORTRAN code and affect running the model are discussed. 
Discussion of the code interfacing the FORTRAN code to UNIX and the X Windowing System, 
written in both FORTRAN and C, will be discussed at a later time. 
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MATHEMATICAL DESCRIPTION 


1. MODEL EQUATIONS 

A. Soil Heat and Moisture Flux Equations 


The slow movement of heat and moisture in a porous medium such as soil can be described by 
diffusion-type equations (Nielson et. al., 1972). In the classical diffusion theory, the flux (the 
amount of a substance crossing a unit area per unit time) is proportional to the negative of the 
gradient of the concentration. The proportionality factor is the diffusion coefficient. The best 
known example of this kind of flow is embodied in Darcy’s Law (Hillel, 1971). The wetness flux 
is 


where q 0 is the flux (cubic centimeters of water per square centimeter per second (cm sec' 1 ), 
K(0) is the hydraulic conductivity (cm sec' 1 ), y(0) is the matric potential (cm), and z is the 
distance from some reference point. The term y - z is the hydraulic head and is the potential 
energy of the soil water (matric plus gravitational energy) per unit weight of water. The function 
y is called the matric potential and is the energy per unit weight required to overcome the 
capillary and adhesion forces that bind the water to the soil. Because work must be done to 
remove water from an unsaturated soil, y is negative. The distance z is the gravitational 
potential energy per unit weight. K and y are functions of volumetric moisture 0 (cm 3 water per 
cm 3 medium). In this application, it is assumed that soil properties change only with depth; thus 
the gradient is a derivative with respect to z. Therefore, Equation 1 may be written as 


% = -m 


dy _1 
dz 




= -K(Q)^l—+K(Q) 
dQ dz 

The second line follows from the chain rule of differentiation. 
Defining a diffusion coefficient 


( 2 ) 


D(8) = K(Q) 


dy 

~dS 


(3) 


yields, when inserted into Equation (2) 

% =-D(G)?L+K(Q) 

dz 


(4) 


The first term in Equation 4 is the diffusion contribution to the moisture flux due to a wetness 
gradient. 
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There is a large body of experimental evidence indicating that thermal gradients induce moisture 
flow (Nielsen et al., 1972). For example, if a uniformly moist soil sample is enclosed in a 
horizontal cylinder and is subjected to a thermal gradient, moisture flows from the warm toward 
the cool end. As field soil temperatures are always changing, an isothermal model such as 
Equation (4) is not complete; a theory that treats both heat and moisture flow in soils is 
necessary. In the following description, diffusion-type expressions for both heat and moisture 
fluxes are presented. The derivation closely follows the work of Philip and De Vries (1957). 
Contributions to heat and wetness fluxes that are proportional to wetness and temperature 
gradients are described. The conservation of mass and energy is the invoked to derive the partial 
differential equations that describe the variation with time of temperature and moisture profiles. 

The diffusive flux of water vapor in a porous medium is modeled as 

Z = (5) 

where q v = vapor flux density (gm cm' 2 sec 1 ) 

D atm = molecular diffusivity of water vapor in air (cm 2 sec' 1 ) 
f = tortuosity and porosity function 
p = density of water vapor (gm cm 3 ) 

Equation (5), with f = 1, describes the diffusion of water vapor in air (Eagleson, 1970). The 
factor f represents the reduced volume available for vapor diffusion in the soil matrix due to 
obstacles such as soil particles and liquid water which adheres to them. An experimentally 
determined graph of f as a function of 0 can be parameterized by a linear function (Jackson et 
al., 1974) 

fiQ) - ct(e-0) (6a) 


where e - soil porosity 

a = constant less than 1 


The diffusivity is a function of temperature and can be adequately modeled by the equation 
(Kimball et al., 1976) 


( 


D =D 

aim o 


273.16 


V 


\ 1.75 


(6b) 


where D 0 = 0.229 (cm 2 sec' 1 ) and T is the absolute temperature. 


The gradient in Equation 5 is to be evaluated in terms of moisture and temperature gradients as 
these are the dependent variables of the model. This can be accomplished by using the 
relationship between vapor density and relative humidity 

P=P.A = P.exp[( V £)/(flr)] < 7a ) 

where p 0 = density of saturated water vapor 
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h = relative humidity 

g = gravitational acceleration constant 

R = gas constant for water vapor = 4.615 X 10 6 (ergs gm' 1 K* 1 ) 

The vapor density p 0 depends on temperature and can be approximated by (Kimball et al. 1976) 

Pi> (T)=cxp[R 0 -(R t rr)] < 7b > 

where Ro = 6.0035 

R, = 4975.9 (K) 

T = temperature (K) 

Equation 7 can be derived from the laws of thermodynamics. Assuming water vapor behaves 
as an ideal gas, an expression can be readily obtained relating the vapor pressure, the 
temperature, and the chemical potential of the gas. The chemical potential and the matric 
potential of liquid water are related because they both represent the free energy of the respective 
phases and the two phases are in equilibrium. The gas density is proportional to the partial 
pressure. 


The gradient in Equation 5 can be expressed in terms of temperature and moisture gradients as 
follows: 


= ^P«/0 = P ,/to + Zrtpp 
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The derivative of h with respect to T can be computed from Equation 7 according to 


( 8 ) 


f \ \dh -\\fg 

T jar ” In* + 


'jYV 


(9) 


-In h 


+ In h 




dy 


The matric potential is dependent on temperature through the dependence of the surface tension 
of water on temperature, which is responsible for the capillary force that binds the water to the 
soil matrix. Therefore, y is proportional to surface tension a (Philip and De Vries, 1957) and 
we may write 


r r 

d\\f _ 



af ” 

v°. 


do 

af 


(10) 


A table giving surface tension at a pressure of one atmosphere as a function of temperature 
(Eagleson 1970) can be fit with the exponential 
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o(T) = a o exp[-7(r-273.16)] 


(ID 


where a 0 = 75.9 (dyne cm 1 ) 
0 = 2.09 x 10 3 (K 1 ) 
T = temperature (K) 


The derivative of \\f with respect to T can be computed using Equations 10 and 11. Equation 9 
therefore is 


dh 

If 


( 


-h \nh 


y + 

v 


The 0 dependence of h is, from Equation 7 


( 12 ) 


dh 

~dQ 


r jJ^L^=h\nh — \n\\f 
RT JdQ dQ 


(13) 


Matric potential \j f typically changes by four to six orders of magnitude over the range of wetness 
values normally found in unsaturated soils. A comparison of Equations (12) and (13) shows that 
the variation of h with 0 is much larger than the variation of h with T, at least over the range 
of temperatures found in soils (273 to 310 K). Therefore, relative humidity h in Equation (8) is 
considered to be only a function of 0. 


Since water vapor behaves approximately like an ideal gas, its density depends primarily on 
pressure and temperature. Therefore, p 0 can be assumed to be a function of temperature only, 
with no dependence on 0. With h depending only on 0, and p 0 depending only on T, Equation 
(8) becomes 





' dh > 

v0 + h 


, 39 , 




Inserting Equation (14) into Equation (5) and using Equation (13) for dh/d0 yields 

c f - ~D t vT - D v0 

* v T yap v t vap 


(14) 


(15a) 


where 


^T. Vap - D a t ^~ Q ) h 


7 ^ 

dT 


J 


(15b) 


D 


D 


Q,vap 


a(E-B)p o gh(dy\ 

d0 


RT 


V 


(15c) 


This is the sought for diffusion expression for the vapor flux. Diffusion coefficients D T vip and 
D 6vip respectively called the ‘’thermal vapor diffusivity*’ and the “wetness vapor diffusivity" 
depend on both 0 and on T. 
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The liquid flux can be computed from "Darcy’s law" (Equation (1)). The gradient of \\r in terms 
of moisture and temperature gradients is 


Y do dr 

Equations (10) and (11) give the derivative of \\f with respect to temperature: 

d\\f _ 


ar 


-YV 


Thus, the liquid flux is 


<7 


(16) 


(17) 

(18a) 


where 


D ar 
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The total moisture flux, q 0 (gm cm' 2 sec 1 ), is the sum of the vapor and liquid fluxes 

This can be written in a diffusion form by adding Equations (15) and (18) to obtain 

<f 0 = -D 0 ?9 - D t *T + K*z 


where 


Dq - liq + D q vap 


Df - Drjiq + ^T.vnp 


(18b) 

(18c) 

(19) 

(20a) 

(20b) 

(20c) 


The volumetric water content, 0, is the volume of moisture per unit volume of soil. Because the 
density of water is 1 gm cm' 3 , 0 also represents the mass of water per unit volume of soil, 
assuming that the water includes the liquid phase plus the gas phase. As 0 represents the mass 
and q e is the mass flux, they arc related by the continuity equation 

£ + **-<> < 21 > 
dt e 

This is a partial differential equation involving 0 as a function of depth and time. An analogous 
diffusion equation can be derived to describe the time dependence of the temperature profile as 
a function of the soil heat flux. 

Fourier’s heat flow equation gives the heat flow due to a temperature gradient as 
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= -**T (22) 

where q^x is the temperature-driven heat flux (calories cm' 2 sec* 1 ) and A is the thermal 
conductivity of the medium (cal cm' 1 sec* 1 K' 1 ) 


To apply this equation to heat transfer in the soil, the effective thermal conductivity of the soil- 
water-air system must be known. A generally accepted model (De Vries, 1975) gives X as a 
weighted average over the thermal conductivity of each soil constituent: 


fX +kfX + k f(X 4J 

t a J a' a vap' (23) 

/ + k f + k f 

J w gJt a J a 

where f w , f, and f t are the volumetric fractions of the liquid, solids, and air constituents 
respectively. (It should be noted that f w and 8 are the same, and the porosity, 8, is equal to f w 
+ f t .) The thermal conductivities of each component are A*, A*, and X,. Factor k, represents the 
ratio of the average thermal gradient in the solid constituents of the soil to the average thermal 
gradient in water. It also depends on the shape and orientation of the soil grains. For spheroid- 
shaped particles, factor k, is given by 


' 3 


1 + 


K 

x. 


- i 


u 


8 t 


K 

X 


- 1 




(24) 


where g s is the shape factor and is equal to 1/2 for cylinders of infinite length, 1/3 for spheres, 
and 0 for disks of infinite radius. The complete model uses a sum over various soil solids, but 
only one representative solid is allowed in this program. The weight factor k, for air can be 
determined from Equation (24) with X equal to the thermal conductivity of dry air. The air shape 
factor g t in this case would have no physical meaning and is usually treated as a variable function 
of moisture that must be determined for each soil type. Therefore the air shape factor k, (and 
not gj is input for the air phase. 


The latent heat absorbed or emitted by the soil as the wetness changes state between the liquid 
and vapor phases can be an important cause of temperature fluctuations. This heat can be 
included in the heat flux by assuming that the vapor flux carries with it a heat flux due to the 
latent heat of vaporization that it absorbed from the soil when it evaporated. This heat flux 
carried by the vapor phase is 

<n,=^v 

where L is the latent heat of vaporization (cal gm* 1 ) and q v is the vapor flux (Equation (15)). 
Both thermal and moisture gradients contribute to q v and therefore contribute to q^ v . The moisture 
contribution is computed by inserting the appropriate term from Equation (15) into the above 
equation to obtain 

ft _ = LD. ?8 (25) 

TA,v(0) Q.vap 

The temperature gradient contribution from Equation (15) is included by increasing the apparent 
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thermal conductivity of the air-filled pores, where the vapor phase primarily exists. This vapor 
has thermal conductivity ^ and carries heat flux -X^VT according to Fourier’s law, where VT 
is the temperature gradient in the pore. However, this heat can also be represented by the 
thermal term in Equation (15) with porosity factor f set equal to 1 as 


-Dh 


f \ 

d l± 

y dT J 


IT 


by equating these two expressions for the same heat flux, the vapor conductivity is found to be 

(26) 


X = LD h 


f \ 
dp 


vap 


K dT ) 


Therefore, the total heat flux in the soil is 


tfk ~ $h,T + ^A,v(0) 

where X is given by Equation (23) and includes the vapor thermal conductivity. 


The total thermal energy per unit volume of medium at temperature T is CT, where C is the 
volumetric heat capacity (cal cm' 3 K 1 ) and T is the temperature in kelvins. The conservation of 
heat energy leads to an equation similar to the conservation of mass for water (Equation (21) 

diCT) + V-X , o (28) 

dt * 

The volumetric heat capacity of the soil is computed as a sum over the capacities of the 
constituents (De Vries, 1975) by 

C=/C+/C+/C (29) 

J s j J w w J a a 


Fractions f s f w f t are the volumetric constituents of solid, water, and air; and C s , Q, and C, are 
the heat capacities of the constituent solids, water, and air respectively. 


Equations (21) and (28) describe the time dependence of soil moisture and temperature profiles. 
In this application, only vertical fluxes are considered; this constitutes a stratified model of the 
soil. Therefore, the gradient operator can be replaced by a derivative with respect to soil depth. 
Thus, the moisture and heat fluxes are, from Equations (20) and (27), given by 


<?e = ~ D i 


?»=-*■ 


f d0 l 

-D t 

r dT ' 



J Z J 


'df* 

K dZj 


- LD , 


6, vap 


-K 
dz 

\ J 


(30a) 


(30b) 


The time derivatives of the moisture and temperature profiles, from Equations (21) and (28), are 
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(31a) 


dS = jK 

dt dz 


dt 


i 


(31b) 


B. Boundary Conditions for the Soil Flux Model Equations 


To solve these equations, boundary conditions must be supplied for both wetness and temperature 
at the air/soil interface and in the bottom layer of the soil profile. In principle, either the fluxes 
q 0 and q* or the variables 8 and T could be specified. In the simulator both heat and moisture 
fluxes q e and are computed at the soil surface. In this way the effects of the environment (ie, 
rainfall, evapotranspiration, radiation, etc.) on the evolution of the soil temperature and moisture 
profiles can be modeled. At the bottom of the profile a mixture of flux and variable boundary 
conditions are used. One can specify constant wetness, a downward wetness flux equal to the 
hydraulic conductivity of the bottom layer, a flux equal to zero, or a constant flux. The 
temperature in the bottom layer is held constant. When both temperature and moisture profiles 
are modeled, the surface fluxes can be found by the solution of the heat balance equation 

5 = R + LE + H (32) 

All fluxes are given as flux densities and are positive downward. S is the heat absorbed by the 
soil, R is the net radiation flux density, L-E is the evapotranspiration heat flux density, and H 
is the sensible heat flux density. This equation can be written with the temperature at the soil 
surface as the only unknown variable. After finding the solution, the soil heat flux density at the 
surface, q^ is set equal to S, and the soil surface moisture flux density q e is set equal to E. 


The heat flux absorbed by the soil can be evaluated by using the discrete form of equation (22) 


5 = -X t 


r T,-T ^ 


where X x = thermal conductivity of the surface soil layer 

Tj = temperature at the center of the surface soil layer 
T, = temperature at the soil surface 
= depth to the center of the first layer 


(33) 


1 . Incident radiation 


The net radiation R is usually divided into average net short and net long wavelength components 
(Eagleson, 1970) 

R = R k + R, (34a) 

short long 
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(34b) 


The net short wave component of R at the surface is 

( 1 -*)(!-*>/« 

where A = surface short wave albedo (constant) 

B = fraction of short wave radiation absorbed by the cloud cover 
\ - insolation at the Earths surface for a cloudless sky 


The insolation at the surface, I^is usually one the the observed meteorological inputs to the 
simulator. However, when MODSOL = 1 the simulator employs an empirical model for I,. 


in which 




( 


7 7 

exp 

n 


< V 


sina 


C 0 + C l 


l°g 10 (since) 
sina 


T 




sina = sin5 sin<}) + cosS cos<(> cost 


(34c) 


(34d) 


where = short wave solar energy flux density incident at the top of the Earth’s atmosphere 
(0.033 cal cm' 2 sec* 1 ) 
n = air turbidity factor (n ~ 2-5) 

a = angle between the sun and the local tangent plane (solar altitude) 

<j> = local latitude 

5 = angle between the sun and the plane of the celestial equator (in radians) 

(-23 deg < 8 > 23 deg) 

x = hour angle of the sun (radians) = W d (t - 12) 
t = hour of the day 

W d = ti/ 12 (rad/hour = 7.2722 x 10 5 rad/sec) 
c 0 = scattering model intercept (0.128) (SCATO) 
c x = scattering model slope (0.054) (SCAT1) 

d = short wave absorption by clouds. During rain storms d = (0.5); else d = CLDABS 
( ATTEN, else CLDABS=CLDIN(I) ) 


The solar declaration, 8, is approximated from the day of the year according to (Bras, 1990). 


8 = 


23.45 7i 
180 


cos 


271 

365 


(172-D) 


(34e) 


In order to rectify the times of simulator output and observations made at a site, simulator noon 
is adjusted to clock noon using an approximation to the equation of time (EQT) (Nautical 
Almanac Office, 1987). Compute first the number of fractional days (FD) into the year in as of 
noon on day D, universal time. 


FD = D + CNOON + [X g (7t/180)(43200/7i)] (24/3600) 


(34f) 



The equation of time, in seconds, for noon at the location of the simulation is approximated by 


EQT = 60 [-7.66 sin((0.9856 FD + TCI ) JL ) 

180 

- 9.78 sin((1.9712FD + TC2) JL)] 

1 80 


(34g) 


and the simulator noon is calculated 


TNOON = CNOON - EQT 

where D = the day of the year (1 < D < 365 or 366) in universal time. This may be the first 
day of the simulation or the day mid way through a multi day simulation; it does 
not change during the simulation. 

FD = fractional days into the year (universal time, in days) 

X g = longitude of the local meridian (west longitude is negative) (in degrees and fractional 
degrees 

TNOON = seconds from integration start to simulator noon 
CNOON = 43200 (number of seconds in 12 hours) 

TCI = first adjustable constant (adjustable for current year) 

TC2 = second adjustable constant (adjustable for current year) 

If the observations are recorded in day light saving time the appropriate adjustment will be made 
by the program. 


When the net incident short wave radiation is to be modeled (when IRNET equals any of 0, 2, 
or 3, or MODSOL=l) the surface albedo, A sur , must be modeled. This is done, when the sun is 
above the horizon, according to a polynomial 

A sur = c t + c 2 sina + c 3 (sina) 2 + c 4 (sina) 3 (34g) 

where the c’s are fit to the surface conditions. Default values are c, = 0.22, c 2 = -0.08, c 3 = c 4 
= 0, which produces an albedo of 0.15 when the sun is 58 degrees above the horizon. When the 
sun is below the horizon, A sur is set to 0. 

The net short wave is = I„ (1 - A, w ). 


The contribution to the net radiation from the long wavelength part of the spectrum is modeled 
by 


- aT> J 


(34h) 
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where T, = air temperature (K) 

e t = the emissivity of the air 
T, = surface temperature (K) 
a = Stefan-Boltzman constant 
e, = emissivity of the surface 

The emissivity of the air is modeled as (Eagleson, 1970) 

e = 0.74 + 0.005 £ (34i) 

a ® 

where e 4 is the vapor pressure of the atmosphere in millibars. The vapor pressure and air 
temperature are supplied by the user, and should be taken from the same height above the 
surface. Equation (34h) is a mathematical statement of the assumption that the air and land 
surface radiate with emissivity of e, and e, respectively. 


It was Camillo’s intention that if a canopy were present, the net short wave radiation (short wave 
absorbed by the canopy) would be computed from an analytical solution to the equation for 
multiple scattering which models short wave radiation in a plant canopy (Camillo, 1987). 
Although this is not yet implemented, the result of Camillo’s development is briefly outlined 
here. 

The single scattering phase function p is the probability of scattering from the incident direction 
(u’,4D to (u,<(>) and is modeled by 

p(w,<J>;u / ,<j> / ) = co [1 + xcos(Q)] 

where co = the single scattering albedo 

x = a measure of ansiotropy in p 
= the scattering angle. 

Direct solar and diffuse sky radiation are integrated seperately over the hemisphere. The 
bidirectional canopy flux is expressed divided by the hemispherically integrated solar radiation 
incoming at the top of the canopy. I is therefore the ratio of the scattered flux to the flux incident 
normal to the surface. 

1(L, M,<t>) =//, 0 (L,m,< i>) + (1 - f)l d (L,u,ty) (35) 

where I(L,u,<t>)= flux of scattered radiation in the canopy at optical depth L in the direction u,<J> 
f = the fraction of diffuse (sky) radiation in the incident short wave flux 
1^= the short wave radiant flux in the canopy from diffuse radiation 
^ = the short wave radiant flux in the canopy from direct solar radiation 
L = the leaf area index (optical depth within the canopy) 
u - the cosine of the zenith angle of the radiation 
<{> = azimuth angle of the radiation 

The directional reflectance due to the diffuse sky radiation flux in (35) is 
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(36a) 


I ko (L,u) = A, ex p 
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and the directional reflectance due to the direct solar flux in (35) is 

l d (L,u,i t>) - A t cxp(k o L) 

+ £ 1 exp(-£ o L)| 

+ K exp 


1 -au 

o 
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1 + k u 


+ h x exp 


-GL 

0 J 

1 

1 + bu 
0 

< U ° J 

( \ 
- GL 

U Gu 
* u ° A 


K U ° > 

U°" 

“* J 


cos(<j) -<J>) 


(36b) 


where G = the average projection of the leaves in the direction u 0 
U 0 , 4> 0 = the direction of the sun 

the constants A 0 , B 0 , A lf Bj, k., a OJ b 0 , h 0 , and \ are all explicit functions of the total 
canopy leaf area index, single scattering albedo, phase function ansiotropy 
parameter, soil albedo, and G(u c ) 


The canopy albedo (the hemispherically integrated reflected radiation over all wavelengths) is 

2k 1 


no) 



(37a) 


and can be evaluated analytically with Equations (35) and (36). Similarly, the flux at the bottom 
of the canopy incident on the soil surface is 


I {L t ) 



(37b) 


where Lf is the total leaf area index (integrated from the top of the canopy to the soil surface). 
The signs indicate radiation direction, positive toward the soil surface. Then, with the soil 
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albedo, A,, (assumed Lambertian), the short wave radiation absorbed by the soil (RS $ ) and canopy 
(RS C ) are 

RS, = (1 - A)(/*(Z. r ) + exp (-GL t /u c ))S 0 (38a) 

RS, = (1 - r(o))S c - RS, ( 38b ) 

where S 0 is the hemispherical integral of the incident flux. The exponential function in Equation 
(38a) accounts for the direct solar radiation which passes through the canopy without being 
absorbed. 

Equations (38a) and (38b) are evaluated for both visible (co-0.2) and near IR (ct>~0.8) wave 
bands. If a two layer model (ICAN-1) is specified, the short wave radiation absorbed by the 
canopy (SLCAN) and soil (SLSOEL) are modelled by the simulator (in SUBROUTINE SUN) 
through the following expressions which must partition 1^, the solar radiation absorbed by the 
surface. 


I soil “ Wt ( "0*039 + 0.13 sina ) (39a) 

U = W ( 0.818 - 0.048 sina ) (39b) 

I Nirr is either observed or computed by the program from the albedo and radiation models in 
SUBROUTINE SUN. These coefficients are current default in the program. It should be noted 
that they depend on the leaf area index of the canopy and must be modified as LAI changes. One 
approach to accomplishing this is through Camillo’s canopy reflectance model (Camilo, 1987). 

The photosynthetically active radiation (PAR) absorbed by the canopy, F, at a particular canopy 
optical depth L, may be estimated by the product of the average projection of leaves in the 
direction u multiplied by the radiation intensity in the downwards direction (negative u), 
integrated over all angles. That is 


2 * 1 



where G(u)= the average projection of the leaves in the zenith direction 8, where u = cos0 

I = Equation (35) evaluated for the value of CD appropriate to photosynthetically active 
radiation (PAR); this is co-0.2. 

With the following approximation for G (Goudriaan, 1977) 

G(u) = <}> 0 + <t>i u 

the integral may be evaluated analytically. In the current program the absorbed PAR is modeled 
(in SUBROUTINE SUN) according to 
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F()p„ = WO. 193 + 0.44 sina) 


( 41 ) 


Again the coefficients must be evaluated for each site according to its leaf area index, perhaps 
with the aid of Camillo’s canopy radiation model. 


The net long wave radiation absorbed by the canopy, RL C , and soil, RL„ are (Ross, 1981) 


RL = (1 -T r )(oeX + ™X - 2oe/ c 4 ) 

RL = T r cseX - a£ X + (1 - 

where T, = the observed air temperature (K) 

T c = the canopy leaf temperature (K) 

T f = the soil surface temperature (K) 

T r = the canopy transmission coefficient for long wave radiation. 


(42a) 

(42b) 


T r is calculated according to 

i 

T = 2 fii du exp(-C(«) L t I u < 43 ) 

0 

where Lj- = the total leaf area index 

G(u)= the average projection of the leaves in the zenith direction 0, where u = cos0 


2. Latent and Sensible Heat Fluxes above the Soil Surface. 

The environment above the soil surface is modelled in two alternative schemes according to 
whether the soil surface is considered to be directly coupled to the atmosphere or coupled in 
series through a plant canopy layer. In the one layer model, the soil and plant canopy surfaces 
are modelled as parallel pathways, both directly coupled to the atmosphere in which the 
meteorological observations are taken. The surface is partitioned into bare soil and canopy 
covered fractions. The sensible and latent heat fluxes from the two fractions are summed 
proportionately. In the two layer model, the soil surface interacts with the plant material of the 
overlying canopy and the canopy air. Canopy surfaces also interact with the canopy air. The 
canopy air, in turn, interacts with the atmosphere above the canopy. Modeling of the radiation 
flux and the resistances to latent and sensible heat flux differs somewhat between the one layer 
and two layer model. In both cases, the surface temperature, T„ which is solved is that 
temperature common to the latent, sensible, and long wave radiation flux models which 
simultaneously satisfies the below ground heat and moisture fluxes and the energy fluxes to the 
atmosphere. 
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a. General models for heat fluxes above the soil surface 
General models for the latent (L-E) and sensible (H) heat fluxes are 


L'E = 

_ P C , e , ~ e a 

(44a) 


Y r, * r. 


H = 

1 1 

Cl 

1 

(44b) 


where the following definitions apply: 

L = the latent heat of vaporization (586 cal gm' 1 of water at 273 K) 

E = the mass of water evaporated (cm/cm 2 ) 

p = air density (1.15 X l(f 3 gm cm' 3 ) 

c p = specific heat of air (0.24 cal gm' 1 ) 

y = psychrometric constant (0.66 mb K 1 ) 

e f = vapor pressure at the surface (mb) 

e, = vapor pressure of the air above the surface (mb) 

T f = temperature of the surface (K) 

T, = temperature of the air above the surface (K) 
r, = aerodynamic resistance (sec cm' 1 ) 
r f = surface resistance (sec cm' 1 ) 

The surface vapor pressure is computed from 

e = h e (45a) 

* sat 

where h is the relative humidity of the surface according to Equation (7a) for a bare soil surface 
or h = 1 for plant surfaces. The saturated vapor pressure at the surface, e SJU , is modelled as a 
function of the surface temperature 

„ _ P.', T . iT ) R t*, T , (45b) 

" 1000 

where R gu is the gas constant for water vapor and is the density of water vapor at saturation, 
Equation (7b). 

b. Expressions for flux resistances 

The variables r, and r, represent the resistance to the diffusion of water vapor from inside the 
evaporating surface to just outside the surface and then into the atmosphere, respectively. 

The aerodynamic resistance to exchange between the surface and the atmospheric boundary layer 
above the surface, r„ may be modelled from the similarity theory of the atmosphere diabatic 
boundary layer wind profile, (see for instance Brutsaert, 1984) as 
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where z = the height of the meteorological data measurements 
z 0 = the roughness length of the surface 
k = von Karman’s constant (0.4) 
d = the zero plane displacement 
Uz - observed wind speed at height z (cm sec' 1 ) 

V F 1 = the stability correction for the wind speed profile 
^2 - the stability correction for the temperature profile 


( 46 ) 


The model may be run without employing a stability correction, implying nearly neutral 
atmospheric conditions (IAERO=0) in which case % and *¥ 2 are set to 0. When I AERO 1 the 
program calculates the stability corrections iteratively using Paulson’s (1970) model. The three 
stability conditions recognized depend on the difference between the air and surface temperatures: 


1. Neutral case, tT f - Tj < 0.1 K, = *P 2 = 0. 


For the other conditions the Monin-Obukhov length, S£, is calculated according to a profile 
formulation 


= 


Tu : 


T 


8<T.-T,) In 


d 

i - — 


( ^ 
1 -y 2 




2. Unstable case, T t - T, < -0.1 K 


x ¥ l = 2 In 


1 +x 


+■ In 


^ 1 + x 2 ^ 


n 


- 2tan _1 (jc) + — 
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where x = (1 - 16 z/S£) 025 

3. Stable case, T a - T, > -0.1 K 
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a) if z M <> 1, then V, = *F 2 = -5(z - zJ/SE 

b) if zM > 1, then % = ¥ 2 = -5(z - zj 


The surface to which the surface resistance r, applies consists of two fractions: bare soil and plant 
canopy surfaces. Separate resistances are modeled for each. 

No generally accepted model of bare soil resistance is known at this time, so the simulator 
employs a function of soil moisture (Camillo and Gurney, 1986). 

r lcil - * (0) * *.<» (<UD - 0(D) ♦ KJP> ■ 1} (47) 

where R,(0), R,(l), and R,(2) = input constants (sec cm' 1 ) 

0 iit (l) = saturation soil moisture content of top soil layer input by the user 
0(1) = the modelled soil moisture in the top soil layer at the time the resistance is 
calculated 

The current default for the soil surface resistance in the program is 0. 


The resistance coefficient for plant canopy surfaces, r^, is modelled as a bulk stomatal resistance 
depending on the photosynthetically active radiation (PAR) absorbed by the canopy and three 
environmental stress factors: vapor pressure deficit, soil moisture deficit, and canopy temperature 
(Sellers 1986, Camillo, unpublished). (See for instance the discussion by Gates 1982). 

Experimental data from cereals and grains support a leaf stomatal conductance (the reciprocal of 
resistance) which is a linear function of the photosynthetically active radiation (PAR) incident 
normal to the leaf (Choudhury and Monteith, 1986) 

_L = — (b. * F) (48) 

r itaf a \ 

where F is the PAR and a t and are empirical constants. The total canopy stomatal conductance 
is the integral of the leaf conductances over the canopy: 

J_ (dL (b l + F(L)) 

a d 

a i 

% (49) 

L r + J F(L) dL 
o 

where the integral may be analytically evaluated (for individual cases). 



19 



The functional forms for the stress factors are not well known. They are included in the following 
way (Charles-Edwards and Ludwig, 1974; Sellers, 1989). 

The stress factor due to vapor pressure deficit is given by 

f(VPD) = 1 - VAPSTR * VPD (50a) 

where VAPSTR = a stress factor input by the user (Default = 0) 

VPD = tJUJ - e mcts 

(Rg. T. p 0 ) / 1000 

e »*t ~ saturation vapor pressure of canopy air 

e mca* = va P° r pressure according to meteorological observation 

p 0 = saturation vapor density from equation 7b using T, 

T, = the air temperature 

The stress factor due to dry soil conditions (high crown potential) is 

W = 1 - expty.* O^cmin - ^c)) (50l>) 

where 'P c = the crown potential computed during the simulation 
^=5 a user selected minimum crown potential 
= a user selected stress coefficient 


The stress factor due to canopy temperature f(T c ) is 


FTMP4 


w(- 10£4 > 

30.0-(-1.0£4) 


X 


1 .0 E4 T cent 
l.0E4-(-\.0E4) 


1TSTEXP 


where 


TSTEXP- 


TSTRH-TSTR0 

TSTRO-TSTRL 


and Tcent * s canopy temperature in Celsius. 


1.0£4-30.0 

30.0-(-1.0£4) 


= 0.994 


(50c) 


These three stress factors are combined 

f(T c ,VPD,'P c ) = f(VPD) fCP c ) f(Tp); f(T c ,VPD, , F c ) > 0.001 (50d) 


The bulk canopy resistance, r^, is then, from Equations (49) and (50c), calculated according to 
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1 


r - 

can 


b x xL T + PARI f(T cy VPD,v c ) 


where 2 l x - 1.2 x 10' 2 (RCANO) 
b t = 6x 10^ (RCAN1) 

PARI = the photosynthetically active radiation absorbed by the canopy (Equation (40)). 


Thus, any increase in one of the three contributors to environmental stress will increase the 
canopy resistance. In the absence of accepted functions for the stress factors, they may all be set 
to zero, f = 1, and rely on fitting the canopy resistance coefficient a,. If a x is not constant from 
day to day at a site, this may reflect need to provide estimated functions for the stress factors. 
(One might take advantage of the expectation that these stress factors vary with soil moisture, 
air temperature, and light conditions to estimate them from available field data.) In the current 
program, fQ = 1. 


c. Surface water contribution to the evaporative flux 

The model allows for water loss by transpiration and by evaporation from the soil surface; 
evaporation from water accumulated on both plant and soil surfaces; and on dew formation. The 
accumulation of surface water from interception and dew formation is expressed as a fraction of 
a maximum possible depth of average accumulation for canopy, x c and soil, x s , respectively, as 

X c = hyh^ (51a) 

X, = hA m « (51b) 


where h c and h, are the respective depths of accumulation. 

These are obtained as solutions to linear, first order differential equations (Equations (52a) and 
(52b)) (Massman, 1983). 

~ = (\-p)R 0 -[D' + d e RJX c + E' c 
at 

^ = pR'-lD' + d'RJX' + E'.-KS 

In Equation (52a) and (52b) the first term is the rate at which the rain is intercepted by the 
corresponding surface, the second term is the rate at which intercepted water drips off the 
canopy, and the third term is the evaporation or dew formation rate. is the rain rate, pis the 
fraction of rain which falls directly through the canopy to the soil surface, and D 0 and are 
empirical constants taken to be 0.12 mm/hr and 0.3 respectively. The fraction p is modelled as 
p = exp(-G(l)/Lr) (Sellers et al., 1986), where G(l) is the average projection of the leaves in the 
direction of the incoming rain. KS is the saturated conductivity of the soil. 


(52a) 

(52b) 
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Loss of standing water from the canopy and soil surfaces is then modelled by 

,, P c , !. ^ew) (53a) 

"V" r X , e<e r (transpiration) 

I ca c i ca 

,,, „ P c , e ." e « . 1. e„e C '(dew) (53b) 

^ X t , e t <e ca (transpiration) 

I ji i 

where r M = the aerodynamic resistance over the canopy 
r s , = the aerodynamic resistance over bare soil 


The surface water heights are updated by 


h t ,c(t + AO = h JC (t) + 


dt 


(54) 


where At is the integrator step size. The heights are constrained to be between zero and their 
maximum allowed values. 


d. The one layer model 

In the one layer model the latent and sensible heat fluxes are computed seperately for the bare 
soil and canopy. The computation for soil and for canopy is each weighted according to the 
fraction of the surface which is bare soil and plant covered. The solved surface temperature, T„ 
common to canopy and soil surfaces, satisfies all the energy budget terms at the soil-air interface. 
The fluxes are then computed using T f . The heat fluxes are given by 

T -T 

H = -pc / — - ( 55a ) 

* r p J * r 


T-T 

H c = 


(55b) 


When there is no standing water on soil or plant surfaces the latent heat fluxes are given by 


L-E = 




/.( 1 -/J 


_Vfi_ 

r ,a + r soU 


(55c) 


pc c, - e 

E'E = -!LZ(l-f t )(\-f) 


r +r . 

ca ecu i 


where f, = the fraction of the surface that is bare soil 

l-f s = the fraction of the surface that is covered by vegetation 
f sw = the fraction of maximum surface water on the soil 


(55d) 
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f cw = the fraction of maximum surface water on the canopy 

l-f sw = the fraction of dry soil surface 

l-f cw = the fraction of dry plant surface. 

e c * = canopy vapor pressure, saturated at canopy temperature 

When there is standing water on soil or plant surfaces the latent heat fluxes are given by 


L'E 

s 


pc 


(55e) 


1 -f) (1 -fj 


where e f * = soil surface vapor pressure, saturated at T, 

e c * = canopy surface vapor pressure, saturated at T, 


(550 


e. The two layer model 


For the two layer model the situation is more complex. The flux network consists of transfers 
between soil surface and canopy air; leaf surfaces and canopy air; and canopy air and the above 
canopy air. The resistance network is comprised, as in the one layer model, of the soil surface 
resistance, r soil , and soil aerodynamic resistance, r M ; bulk canopy resistance, r^, and canopy 
aerodynamic resistance, r M . However, there is the additional aerodynamic resistance between 
canopy air and the outside air, r t . The two aerodynamic resistances r„ and r^ depend on the wind 
speed profile within the canopy. The wind speed and turbulent diffusion coefficient profiles are 
approximated from their values at the top of the canopy, (Choudhry and Monteith, 1988) which 
are estimated according to 


ln 


h-d 


u h=- 


ln 


z-d 


(56a) 




k 2 (h-d) 


ln 


z-d 


(56b) 


23 



where Uh = the wind speed at the top of the canopy 
Uob4 =the observed wind speed 
h = the height of the canopy 
d = the zero plane displacement of the canopy 


z = the height at which the wind speed is observed 
z 0 = the roughness length of the canopy 
k = von Karman’s constant (0.4) 


The aerodynamic resistance between canopy leaf surfaces and canopy air, r^, is 

a a 

^T>j U k 

ca oc 

1 - exp(- “) 

where a = an attenuation coefficient (default = 3.0) 

a = a constant input by user (sec/cm) H (default = 3.0) 

Lr= the (total) canopy leaf area index 

The aerodynamic resistance between soil surface and canopy air, r st , is 


(57) 


r _ h exp(a') 


a'Ki 


exp 




a z 


exp 


f , \ 

-a '(z+d) 


V 


P 


where a’ = a diffusion damping coefficient (= 3.0) 
’ = the roughness length of the soil. 


(58) 


The canopy air temperature, T^, is computed iteratively from equations (60a), (60b), and (60c) 
letting 


H + H, - H = 0 

during the integration interval. Solving for T M gives 


T - 


r r T + r r T+rrT 

ca sa a a sa c a ca i 


r r + r r +rr 

ca ta a sa a sa 


(59) 


A similar procedure using equations (60d), (60e), and (60f) letting L*E C + L-E, - L*E = 0 
permits computing the vapor pressure in the canopy air space. The procedure weights the 
resistances in r, and r 2 according to the fraction of standing water on soil and canopy surfaces. 
The soil resistance r $oil and bulk canopy resistance r^ are obtained according to Equations (47) 
and (50e) respectively. 
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The fluxes into the canopy air space and from the canopy air into the air above are then 
computed according to equations (60a through 60f) 


H = 


pcAT„-T a ) 


(60a) 


// = 


P cAT e - TJ 


(60b) 


H _ = 


P cJJ, - TJ 


(60c) 


L'E = - 


( ^ 

PS 

v* 

( 


(* - o 


/>£ = - 


PS 


£•£ = - 


PS 


|(s - ej 


(s - O 


(60d) 


(60e) 


(600 


where the subscripts c, ca, a, and s refer respectively to the canopy, the canopy air space, the air 
above the canopy, and the soil surface. The subscript r x and r 2 are defined 
r, = r c , for wet soil surfaces 


i i 

= r„ + r soil for dry soil surfaces 
r 2 = for wet canopy surfaces 


r 2 ~ r c« + r cm 


for dry canopy surfaces. 


The latent heat fluxes from wet and dry surfaces are added in proportion to their respective 
fractions in a manner analagous to that used in the one layer model. 


C. Solution of the surface energy budget 

To solve the heat balance equation for T iS one starts by rewriting Equation (33) as (Hillel, 1977) 

T = 5.5 + T, (61) 

' 
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( 62 ) 


Inserting the right hand side of the heat balance equation (32) for S gives 

T = T x + + L-E(T t ) + //(f)} 

The dependence of all three modeled flux terms on the same surface temperature, T„ is thus 
established. This equation, having units of temperature on both sides, is of the form 

T, = F(T.) 

and can be solved by the method of successive approximations. A trial value for T, is chosen, 
FOTJ is evaluated, and a new value for T, results. This procedure can be repeated until 
satisfactory convergence is obtained. In the simulator the air temperature T a is used as a start 
value, and a maximum of four (default) iterations are allowed. Tests have shown that the process 
converges after one or two iterations. Convergence is defined as the absolute value of the change 
in T f between iterations being less than 0. 1 K. The process will always converge if the magnitude 
of the derivative of F is less than 1, 


This derivative is a complicated function of changing meteorological variables, so an analytical 
study of the conditions required for convergence is not feasible. However, from Equation (62) 
it is clear that this derivative is proportional to Zj, the depth to the center of the first soil layer. 


Periods of rainfall can be modelled The user supplies the number of rain storms, the start and 
stop times of the rain (^ and tj), and the total accumulation (r^ for each one. A constant rate 
throughout each time interval is assumed and calculated as 


r - 



(63) 


The short wave radiation attenuation factor during each rain storm must also be supplied. This 
number is equivalent to the cloud attenuation factor B in equation (34b). 


During periods of rain, the evapotranspiration and sensible heat fluxes are set equal to zero (LE 
and H, Equation (32)), and the soil wetness flux at the surface qe, is set equal to the rain rate, 
Equation (63). 


It is possible to remove all temperature dependence from the simulation. (See the description of 
the NAMELIST parameter ITEMPS). In this case the temperature profile, soil heat fluxes, and 
atmospheric heat fluxes are not modelled. However, the evapotranspiration flux must still be 
estimated to provide the surface moisture boundary condition. To do this a Gaussian function of 
time is supplied. 
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(64) 


£(') = -E ma txp[ ~k(l - l m J 7 ] 

where t„„ is the time of maximum demand and E„„ is the rate at this time. The variable k 
which determines the width of the Gaussian can be related to the integrated daily rate as 
follows: 


This gives 


= fdt | E(t) | - £ m „ jdt exp[ -kit - r m J 2 ] 


(65) 


= E 
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k ~ k 
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V 


'day 


( 66 ) 


X / 

The user supplies t^, E m „, and E^y. The simulator computes k from Equation (66), and then 
Equation (65) is used to model the evapotranspiration flux. 


For some simulations it may be simpler to specify the integrated daily total and the width of the 
Gaussian. The maximum rate E^ can be computed from these two. The exponential slope k 
equals 1/t, 2 , where ^ is the time interval between the maximum rate and the time when the rate 
falls to 1/e of this value. Setting Equation (66) equal to 1/t, 2 and solving for E m „ gives 



Therefore the user can compute the required input E mtI from E^y and t e . 

It is also possible to model the surface temperature and the heat balance equation (and thereby 
the effect of the meteorological variables on evapotranspiration) without modeling the soil 
temperature profile. The surface temperature T, and average subsurface temperature T are 
modelled by the force restore method (Lin, 1980). The mathematical model is 

dT s ^ 2S _ 2n_ (T ^ (67a) 

dt a x ' 

= S (67b) 

dt a\J 365 n 

where S is the heat flux absorbed by the soil and 
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a s 

N 71 

In this expression X is the thermal conductivity and c is the heat capacity of the soil surface layer 
and 1 is the number of seconds in a day. The thermal conductivity and heat capacity are 
computed using Equations (23) and (29). 

Since Equation (67) gives the time dependence of the surface temperature explicitly, T, is made 
one element of the state vector and is therefore known. Therefore, no iteration is required to 
solve the heat balance equation. The terms R, LE, and H are evaluated using the state vector 
value for T,, and S is computed from Equation (32). 


D. Root Model 


A model of soil water depletion by plant roots has been included as an extra term in the equation 
of continuity, Equation (31a) 

£. = - *£ - <2(z,0 (68) 

dt dz 

The sink term Q (1/sec) in Equation (68) is positive when water flows from the soil to the plant. 
The mathematical model is (Hillel, 1977; Gardner, 1964) 


Q(z,t) 


<D ,(0,z)-O,(f) 
Q + Q 

s p 


where O p (t) = the crown potential (cm) 

0,(0, z) = the total potential energy of the soil water (cm) 

0,(0, z) = y(0) - z 

where Y = the matric potential 

-z - gravitational potential (z is depth below the surface and is positive). 


(69) 


(70) 


The plant potential O p (cm) in Equation (69) varies with time but is assigned the same value 
throughout the root system. The soil resistance O, (cm-sec) is inversely proportional to the soil 
conductivity and the quantity of active roots 

n = I (7i) 

' B K(Q) P(z) 

where B = constant 

K(0) = hydraulic conductivity (cm/sec) 

P(z) = relative root density at depth z (cm/cm 3 ) 
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The resistance to flow in the roots 


Q p (z) = r/P(z) 


(72) 


where r = specific resistance to flow in the roots (sec/cm). 


Using Equations (71) and (72) for the resistances and Equation (70) for the soil water potential 
energy in Equation (69), after rearranging, gives 


Q - 


BKP[\\f - z - 
1 + BKr 


_ BKP[\f - z - 

Q = Q (73) 

1 + 

Q 

s 

The important model parameters are the relative root density P(z) and the ratio of the resistances 
Qp/H,. No loss of generality results from setting B=l, since its value can be absorbed into the 
definitions of P and r. Since Q is proportional to P, multiplying P at all depths by a constant 
would only change the rate at which the moisture profiles evolve. Since P has the dimensions of 
1/cm 2 , it is commonly thought of as the length of active roots per volume of soil. As yet there 
is no experimental evidence that this is true; the model only requires that P(z) represent the 
relative ability of the roots to absorb water at each depth. The plant potential d> p , commonly 
referred to as the crown potential, is modelled as a response to an atmospheric evapotranspiration 
demand function. 



The discrete model of the sink term as used in this simulator is 


0 - ~ *1 " ^ (74) 

' 1 + r/T 

Q is the value of the sink in the f soil layer, and zj is the depth to the center of this layer. Kj 
and y are the hydraulic conductivity and matric potential of the soil water in the layer. The 
relative root density in each layer Pj and the specific resistance of the plant roots r are input 
parameters. The crown potential T p (t) is modeled as a response to a known transpiration demand 
function E pl (t). The crown potential is computed by requiring that the integral of the sink terms 
over the soil profile be equal to E pl . In its discrete form, this integral is 


*„ = E ga (7S) 

/*! 

where dzj is the thickness of the j* layer and N is the number of layers in the profile model. 
Using Equation (60) for Q and solving for the crown potential gives 
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( 76 ) 


£„(') + - 2)d Zj 

- — - N 

Y, K j p i dz i 

Both B pl and are negative, so T p is also negative. Its magnitude can be large if either the 
demand is large, the soil is dry (so |\j/j| is large), or both. The magnitude of <& p must be less than 
the magnitude of the wilting point O w , which is the largest potential for water intake the plant can 
create before wilting. Thus the crown potential must satisfy the inequality 

d> < d>* < 0 

If <E> p < the simulator will set O p = d> w . Once d> p is calculated the sink term can be evaluated 
for each layer using Equation (74). It must be positive for all layers, to correspond to flow from 
soil to roots. Any of the Qj which are negative are set equal to zero. This procedure is used to 
accommodate experimental evidence that water flow from plant roots to the soil is negligible 
(reference 8). 

The transpiration demand Ep, is computed from equation (60e). 


D. Soil Hydraulic Properties 


Both matric potential and hydraulic conductivity as functions of volumetric wetness may be 
modelled for a wide range of soil types and textures following Clapp and Homberger, 1978 


W) = 


f v* + 3 
_0 


(77a) 


V(0) = V, 


for 


(77b) 


where 0, is the volumetric wetness at saturation, K, and \|/, are the conductivity and matric 
potential respectively at saturation, and b is a parameter determined primarily by the soil texture. 
Representative values are 4 for sand to 1 1 for clay. 


Instead of using this model, the user may choose to provide tables of hydraulic conductivity and 
matric potential as functions of soil moisture. The simulator will perform a linear interpolation 
within the table (or linear extrapolation for soil moisture values outside the range of those 
supplied). An example of such a table is presented in the Appendix. The user specifies the 
choice of model or table look-up via the NAMELIST variable MODHYD. 
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2. METHOD OF SOLUTION (To be added later, unchanged from the original) 


3. PROGRAM SOILSIM: NAMELIST INPUT AND SELECTED OTHER VARIABLES 

(University of Delaware version) 

Variables Definition 
2 May 1990 

Equation numbers refer to the attached SOILSIM Program Description. 

This is the main system input. The NAMELIST name is INPUT, and it is read on unit-5 Subscripts run from 1 
to NL (number of soil layers) unless otherwise indicated. Variables indicated with (*) have been removed from the 
University of Delaware version and those indicated by (**) have been added. 

I. INTEGRATOR CONTROL VARIABLES 


NAME 

TYPE 

DEFAULT 

DESCRIPTION 

TSTOP 

R*8 

8.64D4 

Duration of simulation (seconds). 

TNOON 

R*4 

4.32D4 

(*)Seconds from integration start to noon - times before noon positive, 
afternoon negative. 

HMAX 

R*8 

1.8D3 

Maximum integration step size if IFORCE=l 

IFORCE 

1*4 

1 

Force integration step size to remain less 
than HMAX (0=no, l=yes) 

H 

R*8 

1.0D0 

Initial step size (seconds) if IFORCE=l; 

otherwise H is set to HMAX/512). 

WATERR 

R*4 

1.0E-3 

Moisture error tolerance (E it Eq. 103) 

TEMERR 

R*4 

1.0 

Temperature error tolerance (E i+NL , Eq. 103) 

ED 

R*4 

5.0 

Error window parameter (Eq. 103) 

FI'EMPS 

1*4 

0 

Select soil temperature model: 

0 = no temperature 

1 = model soil temperature profile 

2 = use force restore method 

JBOT 

1*4 

1 

Bottom boundary condition for moisture 

0 = 0 flux 

1 = constant soil moisture in bottom layer 

2 = downward flux equal to hydraulic conductivity. 

3 = constant flux set by user 

QBOT 

R*4 

0.0 

Moisture flux at bottom of profile (cm/sec) if JBOT=3 
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The following variables are for establishing the simulation noon (TNOON) in terms of local mean clock time and 
for calculating the declanation of the sun. 


XLONG 

R*4 

0.0 

(**)Longitude of site in degrees (west longitude is negative) 

XLAT 

R*4 

0.0 

(**)Latitude of the site in degrees (north latitude is positive) 

DAYNUM 

R*4 

0.0 

(**)Day of the year (1 ^ DAYNUM ^ 364 or 365) (Eq. 34e) 

IDST 

1*4 

0 

(**)Clock time is Standard time or Daylight saving time (Eq. 340 
(0 = standard time, 1 = daylight saving time) 

TCI 

1987) 

R*4 

4.02 

(**)First constant in equation of time (degrees) (default values for year 
(Eq. 340 

TC2 

See 

R*4 

17.41 

(**)Second constant in equation of time (degrees) (default for year 1987; 
Almanac for Computers, U. S. Naval Observatory for current 
values) (Eq. 340 


32 



II. OUTPUT CONTROL VARIABLES 


NAME 

TYPE 

DEFAULT 

DESCRIPTION 

DTOUT 

R*4 

1800.0 

Output period (seconds) 

ITABLE 

1*4 

1 

Indicator for amount of print output* 

0 = none 

1 = NAMELIST only 

2 = NAMELIST and boundary conditions 

3 = NAMELIST, boundary, soil conditions. 

NFUNCT 

1*4 

0 

(*)Number of moisture and temperature profiles per plot page. 

NWATRS 

1*4 

0 

(*)Number of soil layers for which moisture is to be plotted as a function 
of time. 

INDXW 

1*4 

10*0 

OI=U0) 

(*)Indices of soil layers for moisture versus time plots 

ITRANS 

1*4 

0 

(*)Plot transpiration flux (0=no, l=yes) 

NWFLUX 

versus 

1*4 

0 

(*)Number of soil boundaries for which moisture flux is to be plotted 
time. 

INDXWF 

1*4 

10*0 

(U=U0) 

(*)Indices of soil boundaries for moisture flux 
versus time plots. 

NTEMPS 

1*4 

0 

(*)Number of soil layers for which temperature is 
to be plotted versus time. 

INDXT 

1*4 

10*0 

(*)Indices of soil layers for temperature versus time plots. 

ISRFT 

1*4 

0 

(*)Plot surface temperature versus time (0 = no, 1 = yes). 

NTFLUX 

1*4 

0 

(*)Index giving number of graphs of heat flux versus time 

ENDXTF 

1*4 

10*0 

(U=l,10) 

(*)Indices of soil layers for which soil heat flux is to be plotted 
versus time. 

IHBAL 

1*4 

0 

(*)Plot components of surface energy balance? (0=no, l=yes) 

IPR 

1*4 

10 

(*)Output unit for printer plots. (Tables are printed on unit 6). 

IUPRT 

1*4 

7 

(**)FORTRAN unit number for formatted output 

BASENAME 

CHARACTER* 1 28 

(**)Pathname and basefile name for input and output files. Extensions 
added to distinguish the files. Input files are ".in", ".met", 
",hyd". Output files are ".out" for meteorological output; ”.prt" 
for formatted output. Additional extensions may be used to 
identify the file type, ie n .ps" for postscript. 
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The following 10 parameters are lower and upper limits for the axes of the printer plots. If a set of lower and upper 
limits are equal, the actual limits will be determined from the data plotted. 


WL, WH 

1*4 

O.O.O.O 

(*)Limits for moisture plots. 

WFL.WFH 

1*4 

0.0, 0.0 

(*)Limits for moisture flux plots. 

TL, TH 

1*4 

O.O.O.O 

(*)Limits for temperature plots. 

TFL, TFH 

1*4 

O.O.O.O 

(*)Limits for heat flux plots. 

HBL, HBH 

1*4 

O.O.O.O 

(*)Limits for heat balance Eq. plots. 

Indices to control accumulation of moisture and heat in soil at boundaries 

NWCUMS 

1*4 

0 

Number of soil boundaries for which 

cumulative soil moisture fluxes are to 
be computed. (<=10). 

IXWCUM 

1*4 

10*0 

(U=l,10) 

Indices of boundaries for cumulative moisture fluxes 

NHCUMS 

1*4 

0 

Number of soil boundaries for which cumulative 
heat fluxes are to be computed (<10) 

IXHCUM 

1*4 

10*0 

(IJ-U0) 

Indices of boundaries for cumulative heat fluxes 

NDISK 

1*4 

0 

Indicator if moisture and temperature profiles 
are to be output to disk (0=no,l=yes). 

IUDISK 

1*4 

12 

Unit number for profile output 


(*) Indicates variable added to the unix version. 
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in. Variables for Defining the Soil Profile 


NL 

1*4 

200 

Number of soil layers (2 <, NL <, 200) 

DZ(I) 

R*4 

200*1.0 

Thickness of soil layers (cm) 

WATER(I) 

R*4 

200*0.25 

Initial volumetric moisture of soil layers. 

TFORCE 

R*4 

2*293 

Initial force restore soil temperatures (Eq. 67) 

TEMPS (I) 

R*4 

200*293.0 

Initial temperature of soil layers (degrees K). 

MODHYD 

1*4 

1 

Source for soil hydraulic functions: 


0 = table look-up 

1 a* Clapp and Homberger model 


IHDUNT 

1*4 

9 

Unit number for look-up tables. 

SATW(I) 

R*4 

200*0.3 

Saturation moisture (0,, Eq. 77) 

SATP(I) 

R*4 

200*-10.0 

Saturation matric potential (cm) (\j/„ Eq. 77b) 

SATK(I) 

R*4 

200*lE-4 

Saturation hydraulic conductivity (cm/sec) (K,, 
Eq. 77a) 

EB(I) 

R*4 

200*5.0 

Soil texture parameter (b, Eq. 77) 

PORSTY(I) 

R*4 

200*0.45 

Soil porosity (e, Eq. 5, 15) 

TCONDS(I) 

R*4 

200*2.5E-3 

Thermal conductivity of soil solids 

(cal cm 1 sec 1 K' 1 ) (X,, Eq. 23) 

VHCAPS 

R*4 

200*0.5 

Volumetric heat capacity of soil solids 
(cal cm' 1 sec' 1 K' 1 ) (C„ Eq. 29) 

SHAPE(I) 

R*4 

200*0.33 

Shape factor (g,, Eq. 24) 

FACTKA 

R*4 

1.4 

Weight factor for air (k,, Eq. 23) 

TCONDW 

R*4 

1.3E-3 

Thermal conductivity of water 

(cal cm' 1 sec' 1 K 1 ) (X*, Eq. 24) 

TCONDA 

R*4 

5.967E-5 

Thermal conductivity of air 

(cal cm' 1 sec 1 K 1 ) (A*, Eq. 24) 

VHCAPW 

R*4 

1.0 

Volumetric heat capacity of water 
(cal cm' 3 K' 1 ) (C w , Eq. 29) 

VHCAPA 

R*4 

3.0E-4 

Volumetric heat capacity of air (cal cm' 3 K 1 ) 
(C„ Eq. 29) 

ALPHA 

R*4 

0.667 

Tortuosity factor (a, Eq. 6a, 15) 
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GNU 


GAMMAO 
GAMMA 1 
IROOTS 
ROOTS(I) 
SPRES 


R*4 

1.0E0 

R*8 

2.09D-3 

R*8 

0.0D0 

1*4 

0 

R*8 

200*0.0D0 

R*8 

1.0D6 


Diffusivity constant (Eq. 15b, 15c) (Looks like a fudge 

factor) 

Surface tension parameter (1/C) (for y in Eq. 11) 

Surface tension parameter for y in Eq. 11) 

Use root model (O=no, l=yes) 

Root density profile (l/cm**2) (P(z), Eq. 71 ,72) 

Root specific resistance (sec/cm) (r, Eq. 72) 
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IV. Variables at the Air/Soil Interface 


Name 

Type 

Default 

Description 

MUNIT 

1*4 

8 

Unit number for meteorological data set 

ICAN 

1*4 


Index for one layer or two layer model 
(0=1 layer, 1=2 layer model) 

(In the two layer model the soil surface 
is coupled to atmosphere serially 
through the canopy air space. Otherwise, 
both canopy and soil surface are coupled 
in parallel to the atmosphere). 

XL AT 

R*4 

45.0 

Latitude (degrees) (4>, Eq. 34d) 

ZHGHT 

R*4 


Height of meteorological data measurements (cm) 
(z, Eq. 46) 

IAERO 

1*4 

0 

Include stability corrections in aerodynamic 
resistance model. 

(0=no, l=yes) 

MAXAER 

1*4 

1 

Max iterations in AERO calculation 

RAERR 

R*4 


Error criterion for iteration in AERO 

DTNEUT 

R*2 

0.2 

Maximum Ts-Ta for neutral conditions in AERO 

NITERS 

1*4 

4 

Maximum number of iterations for HBE solution c 
Ts 

TMPITR 

R*4 

0.1 

Temperature error criterion in HBE solution 

ITRANS 

1*4 

0 

Index to plot transpiration rate 
(0=no, l=yes) 

ML 

1*4 

80 

This variable has not been used to compute 
... looks like page size 

MC 

1*4 

132 

This variable has not been used to compute 
... looks like page size 

THMIN 

R*4 

0.05 

Minimum surface soil moisture which will support 
evaporation. 

SFRAC 

R*8 

0.1D0 

Fraction of ET demand allocated to bare soil evaporation (f, Eq. 55) 

EMAX 

R*8 

3.0D-5 

Maximum rate for evaporation model (cm/sec) 
(E^, Eq. 64, 65, 66). 

EMAXT 

R*8 

4.68D4 

Time of maximum rate (seconds since start of 
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simulation day (t^, Eq. 65). 


EDAY 

R*4 

1.0 

Daily evaporation (cm) (E^, Eq. 65). 

NRAINS 

1*4 

0 

Number of rain storms ( <, 10). 

RNSTRT(I) 

R*4 

10*0.0 

Start times of rain storms (seconds since simulation start) (1=1,10) 
(L, Eq. 63) 

RNSTOP(I) 

R*4 

10*0.0 

Stop times of rain storms (seconds since simulation start) 
(t lt Eq. 63XI=U0) 

RNTOT(I) 

R*4 

10*0.0 

(I=U0) 

Total accumulation for each storm (cm) 
(r^, Eq. 63) 

HSLO 

R*4 

0.0 

Initial soil dew interception (cm) 

HPLO 

R*4 

0.0 

Initial plant dew interception (cm) 

HSLMAX 

R*4 

0.0 

Maximum soil dew interception (cm) (Eq. 51) 

HPLMAX 

R*4 

0.0 

Maximum plant dew interception (cm) (Eq. 51) 

ATTEN(I) 

R*4 

10*0.5 

(I=U0) 

Short wave attenuation factor for each storm (Eq. 34c) 

TURB 

R*4 

2.0 

Turbidity factor (n, Eq. 34c). 

SUNDEC 

R*4 

0.0 

Sun declination (degrees) (8, Eq. 34d) 

EMLONG 

R*4 

0.96 

Long wave emissivity (e„ Eq. 34e). 

RHOVPO 

R*4 

6.0035 

Soil water vapor density coefficient (R„, Eq. 7b) 

RHOVPT 

R*4 

4975.9 

Soil water vapor density coefficient (Rj, Eq. 7b) 

DATMO 

R*4 

0.229 

Diffusion coefficient of water vapor at 0 degrees C 
(cm 2 sec i) (D 0 , Eq. 6b) 

DATMT 

R*4 

0.055E0 

Vapor diffusion model constant (appears unused at present) 

LHEAT 

R*8 

586.0D0 

Latent heat of vaporization of water (cal/gm) 

EMISSO 

R*4 

0.74 

Offset constant for emissivity of atmosphere. (Eq. 34g) 

EMISS1 

R*4 

0.005 

Slope constant for emissivity of atmosphere. (Eq. 34g) 

PSYCHR 

R*4 

0.66 

Psychrometric constant (mb/K) (y, Eq. 53, 55) 

DAIRO 

R*4 

1.29E-3 

Density of air at 0 C used to compute a c p 

DAIR1 

R*4 

0.0045E-3 

Air density model constant (not used in current program) 
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CP 

R*4 

0.24 

Specific heat of air (c pt Eq. 53, 55) 

RSO 

RM 

0.0 

Soil surface resistance offset (sec/cm) (r 0 , Eq 47) 

RSI 

R*4 

0.0 

Soil surface resistance slope (r„ Eq. 47) 

RS2 

R*4 

0.0 

Coefficient of 1/theta in soil surface resistance 
(sec/cm) (Eq. 47) 

VON 

R*4 

0.4 

Von Karman’s constant (k, Eq. 46) 

ZOSOIL 

R*4 

0.2 

Soil roughness length (cm) ( z Eq. 58) 

ZOCAN 

R*4 

25.0 

Canopy roughness length (cm) (z^, Eq. 46, 56, 58) 

DISP 

R*4 

75.0 

Zero plane displacement (cm) (d, Eq. 46, 56, 58) 

V. Parameters to 

i describe the canopy 


RCANO 

R*4 

1.2E-2 

Stomatal resistance constant (a^ Eq. 49, 50e) 

RCAN1 

R*4 

6.0E-4 

Stomatal resistance constant (bj, Eq. 49, 50e) 

RCAO 

R*4 

3.0 

RCA resistance constant ((sec cm' 1 )'*) (a, Eq. 57) 

WNDDMP 

R*4 

3.0 

Wind damping coefficient in canopy model (a, Eq. 57) 

DIFDMP 

R*4 

3.0 

Diffusivity damping coefficient (a’, Eq. 58) 

CPMIN 

R*8 

-1.5D4 

Limiting value of crown potential (cm) Eq. 50b) 

PSISTR 

R*4 

1.0E4 

Crown potential stress coefficient CF^, Eq. 50b) 

VAPSTR 

R*4 

0.0 

Vapor pressure deficit stress coefficient (1/MB) 
(V*. Eq. 50a) 

RSTOMO 

R*4 

0.0 

(not currently used) 

RSTOM1 

R*4 

0.0 

(not currently used) 

TSTRO 

R*4 

30 

Temperature stress coefficient (Eq. 50c) 

TSTRL 

R*4 

-1.0E4 

Temperature stress coefficient (Eq. 50c) 

TSTRH 

R*4 

1.0E4 

Temperature stress coefficient (Eq. 50c) 

CD 

R*4 

0.2 

(not currently used) 

HCAN 

R*4 

100.0 

Canopy height (cm) (h, Eq. 56, 58) 

RLAI 

R*4 

2.0 

Leaf Area Index (L, Eq. 49, 57) 
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EMC AN R*4 1.00 Canopy emissivity (e c , Eq. 42a) 


VI. Parameters for radiation modeling 

IRNET 1*4 3 Index for computing radiation balance, (see METEOR and 

SUN, HBAL, PET) (This needs more work) 

0= no value supplied; program will model incident 
1= net short wave supplied, model computes net long wave. 
2= net long wave radiation supplied 
3= net all-wave radiation supplied, 

4= ???incoming solar supplied, atmospheric long wave modelled 
5= ??? model computes atmospheric long wave. 


MODSOL 

1*4 

0 

Index whether measured or modeled incident solar will be used for 
radiation absorbed by canopy. (0=measured, l=modeled) 

SOIABS 

R*8 

-0.039D0, 

.13D0,2*0.0D0 

Coefficients for solar radiation absorbed by soil 
(Eq. 39a) 

CANABS 

R*8 

0.818D0, Coefficients for solar radiation absorbed by canopy 

-0.048D0, 2*O.ODO (Eq. 39b) 

PARINC 

R*8 

0.193D0, 
0.44D0, 2*0.0 

Coefficients for PAR absorbed by canopy (Eq. 41) 

RLTRA 

R*4 

0.0 

Long wave transmissivity in canopy (T r , Eq.42) 

CALB 

R*8 

.220D0, Canopy albedo coefficients (Eq. 340 

-0.08DO,2*0.0DO 

ALB 

R*4 

0.3 

Surface short wave albedo (Eq.34f) 

CTRANS 

R*4 

0.2 

Cloud transmissivity to short wave (not used in current program) 

CLOUDS 

R*4 

0.0 

Fractional cloud cover (? B, Eq. 34b) (not used in current program) 

SHORTO 

R*4 

0.033333 

Solar constant in incoming short wave model coefficient (cal cm' 2 sec' 1 ) 

SCATO 

R*4 

0.128 

Short wave model constant (Eq. 34c) 

SCAT1 

R*4 

0.054 

Short wave model constant (Eq. 34c) 

PLFRAC 

R*4 

1.0 

Canopy interception fraction ((1-p), Eq. 52) 
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METEOROLOGICAL DATA INPUT 

The meteorological data needed to drive the simulations are input in NAMELIST format The 
NAMELIST name is MET, and the data set is on the unit specified by the parameter MUNIT in the main 
NAMELIST input data set INPUT. 

Each data value (air temperature, vapor pressure, wind speed, and optionally radiation) is used 
throughout a user-specified time interval. Using the parameter IREUSE, one may also specify that the first 24 hours 
of data be reused for succeeding 24 hour periods. 


Following is a description of the data elements in MET: 


NAME 

TYPE 

DEFAULT 

DESCRIPTION 

NMETS 

1*4 

0 

Number of values in each data array (< 1000) 

DELTA 

R*4 

3600.0 

Time interval (seconds) for which each data value is used. 

IREUSE 

1*4 

0 

Reuse first 24 hours of data for succeeding days (0=no, l=yes) 

IROUT 

1*4 

0 

Indicator showing type of radiation data on this data set 


be used by simulator. (Called IRCODE in METEOR and 
IRNET in CANOPY; PET; and SUN). Components not 
supplied with meteorological data will be supplied by models 
(Eq. 34b for net short wave, Eq. 34e for net long wave) 

(0 = none 

1 = Net short wave 

2 = Net long wave 

3 = Net all wave) 

4 = Incoming solar radiation 

5 = ) 


The following are the data arrays, maximum length 1000. 


TMPAIR 

R*4 

1000*293.0 

Air temperature (degrees K) 

VAPORS 

R*4 

1000*15.0 

Vapor pressure (mb) 

WINDS 

R*4 

1000*100.0 

Wind speed (cm/sec) 

SHORT 

R*4 

1000*0.0 

Net short wave radiation (cal cm' 2 sec’ 1 ) 

SLONG 

R*4 

1000*0.0 

Net long wave radiation (cal cm’ 2 sec' 1 ) 

RTOT 

R*4 

1000*0.0 

Net all wave (total) radiation (cal cm’ 2 sec' 1 ) 

CLDIN 

R*4 

1000*0.0 

Fraction of the solar radiation absorbed by the 

cloud cover (B, Eq. 34b). This is used only when solar 
radiation is modelled (IROUT = 0 OR 2) 

SOLAR 

R*4 

1000*0.0 

Downwelling solar radiation (cal cm' 2 sec' 1 ) 
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HYDRAULIC DATA INPUT 


The look-up tables for the hydraulic conductivity and matric potential as functions of the volumetric 
soil moisture are supplied in the NAMELIST data set named HYD. The data set is on the unit number specified by 
the parameter IHDUNT in the main NAMELIST data set, INPUT. 

The interpolation is linear between the nearest moisture values in the table. If the soil moisture 
is either smaller than the smallest table entry or larger then the largest table entry a linear extrapolation is performed. 
To simplify the interpolation logic, the soil moisture values in the look-up table must be evenly spaced. Therefore 


- — 1 ^ 

this option may 

be used only when modelling a homogeneous soil profile. 
The following describes the elements in HYD: 

NAME 

TYPE 

DEFAULT 

DESCRIPTION 

NTHETA 

1*4 

0 

Number of table entries (< 500) 

THETA(I) 

R*4 

none 

Volumetric soil moisture. 

CONDO) 

R*4 

none 

Hydraulic conductivity (cm/sec) 

HE ADO) 

R*4 

none 

Matric potential (cm) 
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DESCRIPTION OF SIMULATOR OUTPUT 


The First two pages of output show values of the variables in NAMELIST INPUT. 

Page 1 has the physical constants which vary with soil depth, and page 2 has values for all other variables. 

OUTPUT AT A PARTICULAR TIME 

TIME elapsed (simulator) time since 0 hours of the first simulation day. 

the format is DDDHHMMSS.MM where 
DDD = days 
HH = hours 
MM = minutes 
SS = seconds 
MM = milliseconds 

CLOCK TIME (**) Local clock time (Standard time or Day Light Saving time) of current output. Hours and 
decimal fractions of an hour. 

CUMULATIVE WETNESS VARIABLES 

TOTAL IN PROFILE is the amout of water (cm) in the profile, calculated from 

NL 

TOTAL = £DZ(i)e(i) 

i-1 


The next lines give integrated moisture fluxes at the boundaries specified by the NAMELIST 
parameters NWCUMS and IXWCUM. The next line shows values of parameters associated with 
the surface energy balance. These are the total, soil, and plant evaporation, and, if IAERO = 1, the 
Monin-Obukhov length, stability corrections % and 4*2, (Equation (46)) bare soil resistance 
(Equation 47), and aerodynamic resistance (Equation 46). The next lines show values of the 
surface energy fluxes. 

SURFACE TEMPERATURE 

If the temperature profile is modelled (ITEMPS=1), then this is the 
solution of Equation 62. If the force-restore method is used (ITEMPS=2) then this is the solution 
of Equation 67a. 

Note that the equations numbers that follow must be updated after Method of Solution is added. 

VALUES OF SOIL VARIABLES 

DEPTH - Depth to the center of the layer 
MOISTURE - Volumetric moisture 

M FLUX - Moisture flux (cm/sec) at the top of the layer (Equation 

78A). The last entry in this column is the moisture flux at the bottom 
of the profile. 

SINK - Q(z,t) (1/sec), (Equation 68) 

HYD COND - Hydraulic conductivity (cm/sec) of moisture in the layer. 

PHEAD - Matric potential (cm) 

DWDT - dG/dt, (Equation 68, 78) 

TEMP - Layer temperature (K) 
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H FLUX 


. Heat flux (cal cm' 2 sec' 1 ) at the top of the layer. The last 
entry is the heat flux at the bottom of the profile. 

TCOND * Layer thermal conductivity (cal cm' 1 sec' 1 K* 1 ). 

VHCAP - Layer volumetric heat capacity (cal cm' 3 ) (Equation 29) 

DTEMPDT - dT/dt (Equation 78c). 


Information at each output time can also be output to disk or tape. The NAMELIST variable NDISK controls 
whether this is done, and IUDISK is the FORTRAN unit number of the DD card which points to the output data 
set. This can be a sequential data set on disk or tape. 

The output records are unformatted. The record length should be at least 8*NL+80 bytes, where NL is the number 
of soil layers. This is 2*NL + 20 words per data record. The records are written by subroutine DSKOUT. They can 
be read with C program convert_outc. 

One header record is output at the start of the simulations, and one data record is output at each simulator output 
time. The header record contains the number of soil layers and thickness of each layer. Each data record contains 
the following: output time, soil moisture in each layer (NL values), temperature in each layer (NL values), surface 
temperature (either the force-restore or that from the heat balance equation solution, depending on whether FILM PS 
is 2 or 1), soil heat flux, net radiation, latent heat flux, sensible heat flux, bare soil evaporation rate, and the plant 
transpiration rate. The record formats are depicted below. 


Structure of Unformatted Output Data Record Created by SOILSIM 
Header Record 


Words Contents 

1 Binary zero 

2 Number of soil layers, Interger*4 

3 Layer thickness, Real*4 

I 

NL+ 2 

NL + 3 Zero, Real*4 

i 

2NL + 2 


Data Record 

Words Contents 
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1-2 


Output time (DDDHHMMSS.SS), RealM 


3 Soil Moisture, RealM 

NL + 2 

NL + 3 Soil Temperature, Real *4 
2NL + 2 

2NL + 3 Surface Temperature (K), TSURF, RealM 

2NL + 4 Net heat + radiation flux at the soil-air interface (cal cm 2 sec l ), SABS, RealM 

2NL + 5 Net short and long radiation summed for soil and canopy (cal cm 2 sec l ), RNET, 

RealM 

2NL + 6 Evapotranspiration (cal cm' 2 sec 1 ), ETHEAT, RealM 

2NL + 7 Sensible heat flux (cal cm' 2 sec' 1 ), SENTOT, RealM 

2NL + 8 Evapotranspiration from soil (cm 3 cm' 2 sec' 1 ), ESOEL, RealM 

2NL + 9 Evapotranspiration from canopy (cm 3 cm' 2 sec' 1 ), EPLANT , RealM 

2NL + 10 Net solar absorbed by soil (cal cm' 2 sec' 1 ), SLSOIL, RealM 

2NL + 1 1 Net solar absorbed by canopy (cal cm' 2 sec 1 ), SLCAN, RealM 

2NL + 12 Net long wave radiation absorbed by soil (cal cm 2 sec 5 ), RNSOIL, RealM 

2NL +13 Net long wave radiation absorbed by canopy (cal cm' 2 sec ! ), RNCAN, RealM 

2NL + 14 Temperature of canopy leaves (K), TCANL, RealM 

2NL + 15 Temperature of canopy air (K), TCANA, RealM 

2NL + 16 Total evapotranspiration from soil and canopy (cm 3 cm' 2 sec' 1 ), ETOT, RealM 

2NL + 17 Sensible heaf flux from soil (cal cm* 2 sec' 1 ), SENSOI, RealM 

2NL +18 Sensible heat flux from canopy (cal cm' 2 sec' 1 ), SENCAN, RealM 

2NL + 19 (**)Modeled net solar radiation absorbed by canopy + soil (cal cm 2 sec ! ), SLTOT, RealM 

2NL + 20 (**)Modeled PAR absorbed by canopy. PARCAN, RealM 
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Other Significant Variables 


RCANOO 

ECANA 

ECANL 

TCANA 

TCANL 

VPDCAN 

SCALE 

VSCALE 

CONI 

CON2 

HSL 

HPL 

SLRAT 

SLRAT1 

PLRAT 

PLRAT1 

CPI 

CP2 

VAPOS 

VAPSOI 

VAPZO 

DA VP 

ESODL 

SENSOI 

SENCAN 


Canopy stomatal resistance (r rt , Eq. 49) 

Canopy air vapor pressure 

Canopy air space vapor pressure 

Canopy air temperature 

Canopy leaf temperature 

Canopy air vapor pressure deficit 

Ratio of plant potential transpiration to actual transpiration 

Scale factor of 1/1 000 used in Eq. 45b. 

pCp/lLy) 

PS 

Interception on soil surface (Eq. 52b) 

Interception on plant canopy surfaces (Eq. 52a) 

Ratio of depth of surface water to its maximum depth on soil surface (Eq. 51b). 
Complement of LSRAT (1 -SLRAT) 

Ratio of depth (HPL) of water on plant surfaces to its maximum depth (HPLMAX) (Eq. 
51a. 

Complement of PL ART (1 -PLRAT) 

CONI * plant cover fraction * SLRAT / RCA (Eq. 55) 

CONI * plant cover fraction * SLRAT1 / (RCA + RCAN) (Ea. 55) 

Saturation vapor pressure of soil surface layer 
Actual calculated vapor pressure of the soil 
Observed vapor pressure of atmosphere 
The difference VAPSOI - VAPZO 
Evaporation from soil 
Sensible heat flux from soil 
Sensible heat flux from canopy 
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CPOT 

TCENT 

TO 

TSTEXP 

SINA 

CAN ALB 

CANM1 

RSOLAR 

RSHORT 

RLONG 

RNET 

RNCAN 

RNSOIL 


Crown Potential, (Eq. 69) 

Temperature Celsius 
273.16 K, 0° Celsius 

Exponent on canopy stress temperature function 
Sine of the solar altitude angle 
Canopy albedo (CALB(l) + SINA) 

Canopy absorption coefficient (1 - CANALB) 

Downwelling solar radiation 

Net short wave radiation 

Net long wave radiation above canopy 

Net allwave radiation above canopy 

Net long wave radiation from canopy 

Net long wave radiation from soil surface 
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